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ABSTRACT 

In  order  to  determine  the  best  approach  for  the  modeling  of 
tracked  vehicle-terrain  interactions,  a  comparison  is  made 
between  terramechanics  and  continuum  mechanics 
plasticity  soil  models  that  can  be  used  in  vehicle  dynamic 
simulations.  The  absolute  nodal  coordinate  formulation 
(ANCF)  which  can  be  used  in  multibody  system  (MBS) 
dynamics  to  model  large  rotations  and  large  deformations 
will  provide  a  novel  framework  for  this  comparison.  First,  a 
brief  review  of  the  analytical  derivation  and 
implementation  of  a  terramechanics  soil  model  for  drawbar 
pull-slip  and  pressure-sinkage  are  presented  in  order  to 
establish  the  basic  assumptions  underlying  this  model.  An 
assumed  shape  and  mass  for  the  vehicle-terrain  interface 
and  an  empirically  derived  constitutive  model  are  two  clear 
assumptions  inherent  in  terramechanics  based  approaches. 
The  (modified)  Bekker  soil  models  are  integrated  finite 
element  (FE)  ANCF  interpolations  to  determine  the 
generalized  forces.  Continuum  mechanics  soil  models  that 
are  suited  for  tracked  vehicle-terrain  interaction  are 
identified  and  integrated  into  the  internal  force  calculation 
of  ANCF  MBS  algorithms.  The  paper  discusses  important 
fundamental  issues  that  must  be  addressed  when 
implementing  continuum  mechanics -based  soil  models  as 
well  as  terramechanics  models  into  MBS  algorithms  for 


modeling  complex  tracked  vehicle/soil  interactions.  The 
improvement  of  the  vehicle-terrain  interface  estimation 
resulting  from  FE  methods  can  provide  a  terramechanics 
approach  which  may  be  scalable  to  differing  vehicles  and 
loadings.  It  is  found  that  the  assumed  existence,  uniqueness, 
and  evolution  of  the  yield  surface(s)  of  continuum  soil 
models  play  a  critical  role  in  predicting  the  soil  behavior 
while  also  providing  a  rational  method  for  improvement  of 
soil  constitutive  modeling. 

1  INTRODUCTION 

In  the  study  of  vehicle-terrain  interactions,  there  are 
two  common  methods  and  models  used:  terramechanics 
and  continuum  mechanics  soil  models.  Broadly  speaking, 
terramechanics  is  the  study  of  the  relationships  between  a 
vehicle  and  its  environment.  Some  of  the  principal 
concerns  in  terramechanics  are  developing  functional 
relationships  between  the  design  parameters  of  a  vehicle 
and  its  performance  with  respect  to  its  environment, 
establishing  appropriate  soil  parameters,  and  promoting 
rational  principles  which  can  be  used  in  the  design,  and 
evaluation  of  vehicles  (Wong,  2010).  The  continuum 
mechanics  approach  is  comparatively  new  in  the 
application  to  the  study  of  vehicle-terrain  interaction.  This 
method  is  of  special  advantage,  since  the  deformations  of 
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the  wheel,  tire  and  ground  can  be  captured  as  the  result  of 
their  mutual  interaction.  The  continuum  mechanics 
approach  is  especially  suitable  for  clearing  up  the 
phenomena  concerning  deformations,  such  as  the  effect  of 
inflation  pressure,  or  bulldozing  resistance.  Additionally, 
the  continuum  mechanics  approach  provides  new 
possibilities  to  investigate  the  pressure  transfer  onto  and 
into  the  ground,  and  therefore  on  soil  compaction  (Schmid, 
1995).  Terramechanics  models  can  be  categorized  into 
empirical  and  analytical  models. 

Empirical  terramechanics  models  use  established 
experimental  measurements  of  appropriate  parameters, 
properties,  and  behaviors  of  soil;  these  experimental  results 
are  then  used  to  establish  empirical  relationships  that  could 
be  used  to  predict  at  least  qualitatively  the  response  of  soils 
under  various  conditions  (Bekker,  1969).  Parametric 
models,  which  are  based  on  experimental  work  and  have 
been  widely  used,  offer  practical  means  by  which  an 
engineer  can  qualitatively  evaluate  tracked  vehicle 
performance  and  design. 

Research  in  the  field  of  terramechanics  has 
continued  over  many  decades.  Since  the  tractive 
performance  of  a  vehicle  depends  significantly  on  the 
normal  and  shear  stress  distribution  at  the  soil-vehicle 
interface,  many  attempts  have  been  made  to  predict  such 
distributions  beneath  tracks  and  wheels  (Schmid,  1995). 
The  following  gives  a  brief  overview  of  the  developments 
in  terramechanics  and  continuum-based  soil  models. 

Bekker  (1956,  1960)  first  proposed  a  theoretical 
study  regarding  the  pressure-sinkage  relationship  of  terrain, 
as  well  as  the  response  to  repeated  loading  and  unloading. 
Janosi  and  Hanamoto  (1961)  proposed  a  model  for 
predicting  the  ground  pressure  distribution  beneath  a  rigid 
track.  They  assumed  the  normal  stress  distribution  to  be 
higher  at  the  rear  end  of  the  track.  The  predicted  ground 
pressure  distribution  was  trapezoidal  with  higher  pressure 
at  the  rear  end  of  the  track.  Reece  (1965)  improved 
Bekker's  model  by  making  the  parameters  dimensionless. 
This  single  equation  could  then  account  for  different  plate 
shapes. 

Garber  and  Wong  (1981)  developed  analytical 
methods  for  predicting  the  pressure  distribution  beneath 
tracks  taking  into  account  all  major  design  parameters  of 
the  vehicle  as  well  as  the  pressure-sinkage  characteristics  of 
the  soil.  Wong  et  al.  (1984)  continued  this  work  and 
established  a  model  for  the  prediction  of  the  distribution  of 
ground  pressure  and  the  tractive  performance  of  a  tracked 
vehicle,  including  the  effects  of  repetitive  loading.  Wong 
and  Huang  (2006)  evaluated  the  tractive  performance  of 
wheeled  and  tracked  vehicles  assuming  that  the  tracks  on  a 
tracked  vehicle  have  an  essentially  flat  and  rectangular 
contact  area  with  uniform  normal  pressure  and  the  same 
contact  length,  as  well  as  assuming  the  vehicle  weight  is 


uniformly  distributed  among  the  tracks. 

Similarly,  Okello  (1994,  1998)  suggested  a 
traction  model  that  takes  into  consideration  soil  recovery 
for  the  repetitive  loading  occurring  under  road  wheels  and 
rubber  tracks.  Rubinstein  (2007)  modified  Bekker’s 
original  model  by  adding  a  viscous  friction  element  which 
allows  modeling  the  dependence  of  the  normal  pressure  on 
the  sinkage  velocity.  Recently  Irani  (2011)  has  proposed  a 
new  form  for  the  pressure  sinkage  relationship  based  on  the 
expression  proposed  by  Reece  (1965)  to  capture  the 
dynamic  oscillations  observed  for  a  wheel  with  grousers. 

Similarly,  research  into  continuum-based  soil 
models  has  continued  over  many  decades.  At  present,  there 
are  many  continuum  mechanics-based  soil  models  that 
employ  different  assumptions,  of  which  some  are  suited  for 
finite  element  implementations.  The  following  gives  an 
introduction  to  some  of  the  more  common  continuum- 
based  single-phase  plasticity  models. 

The  Mohr-Coulomb  model  is  one  of  the  oldest 
and  best-known  models  for  an  isotropic  soil  (Goldscheider, 
1984).  Initially  the  yield  surface  was  used  as  a  failure 
envelope,  and  still  is  in  geotechnical  practice.  It  was  later 
adopted  as  a  yield  surface  for  plasticity  models.  In  two 
dimensions,  the  yield  surface  of  the  Mohr-Coulomb  model 
is  defined  by  a  linear  relationship  between  shear  stress  and 
normal  stress.  The  failure  envelope  defined  by  the  Mohr- 
Coulomb  model  includes  discontinuous  slopes  between 
failure  surfaces. 

A  simpler  method  to  handle  the  discontinuities  is 
to  use  a  smooth  approximation  to  the  yield  surface.  Drucker 
and  Prager  (1952)  initially  proposed  a  cone  in  principal 
stress  space  ,  by  adding  a  pressure-dependent  term  to  the 
classical  von  Mises  yield  surface.  Like  von  Mises  plasticity, 
one-step  return-mapping  can  be  achieved  for  linear 
hardening,  making  the  model  quite  efficient  to  implement. 
While  the  associative  model  over  predicts  dilatation,  non- 
associative  versions  correct  this  (Drucker  et  al.,  1952). 
Initially  developed  as  an  elastic-perfectly  plastic  model,  i.e. 
with  no  change  in  the  yield  surface  on  loading,  researchers 
later  added  hardening  of  the  yield  surface  parameters  to  the 
model  in  various  forms.  See,  for  example,  Vermeer  and  de 
Borst  (1984)  for  a  relatively  sophisticated 
phenomenological  hardening  model. 

The  modified  Cam-Clay  (MCC)  model  by 
Roscoe  et  al.  (1963)  is  based  on  critical  state  theory  and 
was  meant  to  capture  the  properties  of  near-normally 
consolidated  clays  under  triaxial  compression  test 
conditions.  The  yield  surface  is  assumed  to  have  an 
elliptical  shape  that  may  be  expanded  with  the  increase  of 
volumetric  strain.  Cam-Clay  models  can  predict  failure  and 
the  nonlinear  stress-path  dependent  behaviors  prior  to 
failure  accurately,  especially  for  clay  type  soils  (DiMaggio 
and  Sandler,  1971). 
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Cap  plasticity  models  were  developed  to  address 
the  shortcomings  of  the  Cam-Clay  type  models,  especially 
for  high-pressure  behavior.  Drucker  et  al.  (1957)  first 
proposed  that  “successive  yield  surfaces  might  resemble  an 
extended  Drucker-Prager  cone  with  convex  end  spherical 
caps”  (Chen  and  Baladi,  1985).  As  the  soil  undergoes 
hardening,  both  the  cone  and  the  end  cap  expand.  This  has 
been  the  foundation  for  numerous  soil  models. 

The  Sandia  GeoModel  builds  on  the  Cap  model 
with  some  modifications.  It  is  capable  of  capturing  a  wide 
variety  of  linear  and  nonlinear  model  features  including 
Mohr-Coulomb  and  Drucker-Prager  plasticity  depending  on 
the  model  parameters  incorporated.  Unlike  the  Cap  model, 
the  cap  surface  and  shear  yield  surface  are  connected  in  a 
smooth  manner,  and  the  model  also  accounts  for 
differences  in  triaxial  extension  and  compression  strength 
using  either  a  Gudehus  or  William-Warnke  modifying 
function. 

Plasticity  models  such  as  those  described  above 
do  not  include  strain-rate  dependent  behavior  often 
observed  in  soils  under  rapid  loading.  These  viscous  effects 
are  more  pronounced  in  the  plastic  region  of  most  clay  soils 
and  rate-independent  elastic  response  is  generally  adequate 
for  practical  engineering  applications  (Perzyna,  1966; 
Lorefice,  2008).  The  models  described  above  can  be 
modified  to  account  for  rate-dependent  plastic  effects.  Such 
viscoplastic  models  are  more  accurate  under  fast  loading 
conditions.  However,  it  is  difficult  to  determine  the  correct 
value  of  the  material  time  parameter  if  the  stress  history  is 
not  known. 

Soil  is  not  always  an  isotropic  material.  Layering  and 
fracture  networks,  as  well  as  compaction  and  other  history 
effects  may  give  the  soil  higher  strength  or  stiffness  in 
certain  directions.  Often  the  effects  impart  different 
strength  and  stiffness  in  one  plane.  Anisotropy  may  also  be 
addressed  using  fabric  tensors  (Wan  and  Guo,  2001).  Other 
anisotropic  models  include  the  work  of  Whittle  et  al.  (1994), 
and  the  S-CLAY  1  model  (Wheeler  et.  al.,  2003),  which 
builds  on  the  MCC  model.  Detailed  review  of  anisotropic 
soil  models  is  beyond  the  scope  of  this  article,  however, 
and  the  reader  is  referred  to  the  above  references. 

This  paper  aims  to  compare  an  existing 
terramechanics  soil  model  with  a  continuum  mechanics- 
based  soil  model.  A  brief  review  of  the  analytical  derivation 
and  implementation  of  a  Bekker's  terramechanics  soil 
model  for  pressure-sinkage  and  drawbar  pull-slip  is 
presented.  Their  suitability  for  incorporation  into  FE/MBS 
simulation  algorithms  is  also  discussed.  Section  2  outlines 
the  analytical  approaches  used  in  terramechanics.  Section  3 
describes  the  continuum  mechanics-based  soil  models. 
Section  4  describes  the  ANCF  implementation.  The 
comparison  of  the  two  models  is  given  in  Section  5. 
Section  6  offers  a  summary  and  describes  the  direction  of 


future  work. 


2  TERRAMECHANICS  BASED  SOIL 
MODELS 

Three  common  formulations  of  the  pressure-sinkage 
relationship  are  widely  used  in  vehicle-terrain  interaction 
investigations. 

The  first  common  pressure-sinkage  relationship 
for  mineral  terrains  was  proposed  by  Bekker  (1960) 


P  =  {kJb+k<p)zn  =keqzn  (1) 

where  p  is  pressure,  b  is  the  radius  of  a  circular  plate  or  the 
smaller  dimension  of  a  rectangular  plate,  n  is  a  non- 
dimensional  soil  constant  related  to  the  characteristics  of 
the  soil  being  modeled,  kc  and  ktjl  are  pressure-sinkage 
parameters  related  to  the  cohesion  and  the  angle  of  shearing 
of  the  material,  and  z  is  the  sinkage. 

Besides  what  has  been  introduced.  Equation  (1)  is 
now  widely  used  in  track-terrain  interaction  studies.  For 
examples  of  such  studies  refer  to  Wills  (1963),  Ryu  et  al. 
(2003),  Garber  (1984),  Okello  (1994,  1998),  Rubinstein  et 
al.  (2007),  Park  et  al.  (2008).  Similarly  for  tire-terrain 
interactions  see  Mao  (2008),  Sandu  et  al.  (2010), 
Schwanghart  (1991),  Harnisch  (2005),  Gibbesch  and 
Schafer  (2005),  Li  and  Sandu  (2007),  Meirion-Griffith  and 
Spenko  (2011)  ,  Li  and  Sandu  (2007),  and  Ishigami  et  al. 
(2007). 

Another  common  formulation  of  the  pressure- 
sinkage  relationship  was  proposed  by  Reece  (1965) 

P  =  {ck'c  +  (2) 

where,  kc  and  ki  are  dimensionless  soil  parameters  which 
replace  Bekker's  kc  and  k(j  parameters.  Equation  (2)  can 
be  used  on  wheel-terrain  interaction  such  as  in  Irani  (2011) 
and  tire-terrain  interaction  such  as  in  Senatore  (2011). 

Similarly,  among  the  more  common  formulations  of 
the  pressure-sinkage  relationships  are  (Reece,  1965) 


-  cos  0i 


(4) 


where  cr,  is  the  forward  part  of  the  normal  stress,  a,  is 
the  rear  part  of  the  normal  stress,  n  is  the  sinkage  exponent, 
h  and  ki  are  pressure  sinkage  moduli,  b  is  the  wheel  width, 
r  is  the  wheel  radius,  0X ,  0 m  are  the  entry  angle  and  the 
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angular  location  of  the  maximum  normal  stress, 
respectively.  Equation  (3)  and  (4)  can  be  used  in  wheel- 
terrain  interactions  such  as  in  Wong  et  al.  (1967),  Ding  et  al. 
(2009),  and  Shibly  et  al.  (2005).  For  a  flexible  tire  acting  on 
deformable  terrain.  Equation  (3)  and  (4)  can  also  be  used  if 
radius  r  is  taken  to  be  the  effective  radius  (Senatore  et  al., 
2011). 

For  'plastic'  soils  which  do  not  exhibit  a  'hump'  of 
maximum  shear  stress,  the  following  modified  version  of 
Bekker's  equation  containing  only  one  constant  was 
proposed  by  Janosi  and  Hanamoto  (1961)  and  is  widely 
used  in  practice: 


J/K  (5) 

where  K  is  usually  referred  to  as  the  shear  deformation 
parameter.  It  is  a  measure  of  the  magnitude  of  the  shear 
displacement  required  for  the  development  of  maximum 
shear  stress. 

It  has  been  found  that  the  Mohr-Coulomb  equation 
describes  the  normal  stress-shear  stress  displacement 
relation  adequately  in  many  cases 

•W  =c  +  ptm<f>  (6) 

In  this  equation,  s  is  the  maximum  shear  stress,  p 
is  the  normal  stress,  and  c  and  (f)  are  the  cohesion  (or 
adhesion)  and  angle  of  shearing  resistance,  respectively.  In 
deriving  the  values  of  c  and  (f)  from  measured  shear  data, 
the  effect  of  grouser  height  and  of  the  shear  ring  should  be 
taken  into  consideration  (Reece,  1964). 

The  tractive  performance  of  a  link  can  be 
predicted  once  the  normal  pressure  and  shear  stress 
distribution  are  known  at  a  given  slip,  it  is  usually 
characterized  by  its  motion  resistance,  tractive  effort,  and 
drawbar  pull  (the  difference  between  tractive  effort  and 
motion  resistance)  as  functions  of  slip. 

If  it  is  assumed  that  a  track  link  has  a  rectangular 
contact  area  with  uniform  normal  pressure,  the  same 
contact  length  /,  and  the  link  and  terrain  is  horizontal,  then 
under  steady-state  operating  conditions,  the  link  tractive 
effort  F  and  drawbar  pull  Fd  at  a  given  slip  i  can  be 
given  by  (Bekker,  1969) 


of  track  link  in  contact  with  the  terrain,  p  is  the  normal 
pressure,  and  R  is  the  motion  resistance  of  track  link. 


3  CONTINUUM  MECHANICS-BASED  SOIU 
MODELS 

The  original  Cam-Clay  model  has  not  been  as  widely 
used  for  numerical  predictions  as  the  modified  Cam-Clay 
(MCC).  The  qualifier  “modified”  is  often  dropped  when 
referring  to  the  modified  Cam-Clay  model  (Wood,  1990). 
Cam-Clay  models  can  predict  failure  and  the  nonlinear 
stress-path  dependent  behaviors  prior  to  failure  accurately, 
especially  for  clay  type  soils  (DiMaggio  and  Sandler,  1971). 
The  function  for  the  yield  surface  of  the  infinitesimal  MCC 
model  is  defined  as 


q2-M2 


=  0 


(10) 


Here,  q  is  the  deviatoric  stress,  p  is  the  effective  mean 
stress,  the  pre -consolidation  stress  pc  acts  as  a  hardening 
parameter,  and  the  stress  ratio  M  —ql  p  at  critical  state  is 
related  to  the  angle  of  friction  through  the  relationship 
M  =  6sin(0)/(3  +  sin(0))  .  Cam-Clay  models  began  with 
infinitesimal  strain  assumptions  and  have  been  developed 
to  the  case  of  finite  strains  in  Boija  and  Tamagnini(1996). 
Furthermore,  Cam-Clay  models  have  been  extended  to 
capture  the  cyclic  behavior  of  soils  (Carter  et  al.,  1982). 

The  constitutive  equations  that  govern  the 
behavior  of  the  hyperelastic  elastoplastic  finite  deformation 
Cam-Clay  model  (Borja  et  al.,  1996)  are  summarized  as 
follows.  The  yield  function  is  defined  by 

f  =  f{P,Q,Pc)  =  ^  +  P{P-Pc)  (11) 

where  P,  Q,  Pc,  and  M  are  the  finite  deformation  analogs  of 
the  parameters  defined  for  the  infinitesimal  case  using  the 
Kirchhoff  stress  tensor.  The  hardening  law  expressed  in 
terms  of  the  plastic  component  of  the  volumetric  strain  is 
given  by 


F  =  bj  (c+/?tan^)(l-e“/K)<ir 
=  £>/(c+/?tan^)  1  -  (1  -  e~'l/K ) 


(8) 


Fd=F-R  (9) 


where  b  is  the  contact  width  of  the  track  link,  l  is  the  length 


4 

P 


=  -©C 


jp 

Jp ’ 


0  =  - 


1 


A  —  k 


(12) 


where  J p  is  the  plastic  component  of  the  Jacobian.  The 
parameters  A  and  k  can  be  calculated  from  the 
corresponding  infinitesimal  model  analogs.  The  discrete 
flow  rule  at  time  tn+1  for  implicit  time  integration  in  the 
space  defined  by  the  elastic  Eulerian  logarithmic  stretches 
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can  be  written  as 


=  £ 


e  t  r 
n+ 1  “ 


dg_ 

n+bp 


(13) 


where  /?  is  the  Kirchhoff  stress  tensor,  and  A</>  is  a 
plastic  multiplier. 

The  above  equations  can  be  shown  to  lead  to  the  following 
set  of  equations  that  can  be  used  to  define  a  scalar  return 
mapping  algorithm  (Borja,  1998)  in  the  invariants  of  the 
elastic  logarithmic  stretches 


f{P,QPc)<  0  ,A</>  > 


A  <Pf(PQP)  = 


(14) 


significant  change  in  geometry  may  require  the  use  of 
elements  that  employ  gradients  and  accurately  capture 
curvature  changes.  This  requirement  can  be  met  using  the 
FE  absolute  nodal  coordinate  formulation  (ANCF). 

ANCF  finite  elements  do  not  employ 
infinitesimal  or  finite  rotations  as  nodal  coordinates; 
instead,  absolute  slopes  and  displacements  at  the  nodal 
points  are  used  as  the  element  nodal  coordinates.  The 
position  vector  r  7  of  an  arbitrary  point  on  element  j  can 
be  defined  in  a  global  coordinate  system  XYZ  as 
r7  =  S7  (v7  ,y',  z7e7  (f).  In  this  equation,  x  ' ,  y' ,  and 
ZJ  are  the  element  spatial  coordinates,  S7  is  the  shape 
function  matrix,  e7  is  the  vector  of  element  nodal 
coordinates,  and  t  is  time.  The  nodal  coordinate  vector  e  J,‘: 
at  node  k  can  be  defined  as  follows  (Shabana,  2005) 


An  example  implicit  integration  scheme  for  the  finite 
deformation  Cam-Clay  plasticity  soil  model  can  be 
developed  by  considering  a  set  of  simultaneous  nonlinear 
equations.  An  application  of  the  Newton-Raphson  method 
can  be  used  to  solve  this  set  of  nonlinear  equations.  To  this 
end,  the  residual  vector  r  and  the  vector  of  unknowns  x  are 
written  as  follows: 


r 


p 


'+rAt/>  — 

v  v  dP 

dQ 


f 


(15) 


The  Newton-Raphson  solution  procedure  requires  the 
iterative  solution  of  the  algebraic  system 
\drp  /dxppsx.p  =-rp  ,  where  Ax''  is  the  vector  of 
Newton  differences.  A  closed  form  expression  for  the 
consistent  tangent  operator  \drp/dxp)  can  be  found  and 
the  algorithm  can  be  made  more  efficient  by  the  application 
of  the  static  condensation  technique  (Borja  et  al.,  1996). 


4  ANCF  IMPLEMENTATION 

The  FE  implementation  of  the  soil  mechanics 
plasticity  equations  requires  the  use  of  an  approach  that 
allows  employing  general  constitutive  models.  The 
vehicle/soil  interaction  can  lead  to  a  significant  change  in 
geometry  that  cannot  be  captured  using  finite  elements  that 
employ  only  translational  displacement  coordinates  without 
significant  refinement.  In  some  soil  applications,  such  a 


Fully  parameterized  ANCF  finite  elements  allow 
using  a  general  continuum  mechanics  approach  to  define 
the  Green-Lagrange  strain  tensor  8  =  (j  T  J  - 1 J/2  ,  where 
J  is  the  matrix  of  position  vector  gradients.  In  dynamic 
soil  problems,  ANCF  leads  to  a  constant  inertia  matrix  and 
to  zero  Coriolis  and  centrifugal  forces.  The  mass  matrix 
obtained  using  ANCF  finite  elements  can  always  be  written 
a  Mj  =  l pJSJ  SJdVj  ,  where  p7  and  V7  are, 
respectively,  the  initial  mass  density  and  initial  volume  of 
the  finite  element.  ANCF  finite  elements  allow  for  straight 
forward  implementation  of  general  constitutive  models 
including  the  continuum  mechanics-based  soil  model 
discussed  in  this  paper. 

For  a  finite  element  or  deformable  body,  the  principle 
of  virtual  work  can  be  written  using  the  reference 
configuration  as 

J  prTSrdV  +  Jo,,, :  6zdV  -  J  f^SrdV  =0  (17) 

V  V 

In  this  equation,  V  is  the  initial  or  reference  volume,  p 
is  the  initial  mass  density,  r  is  the  global  position  vector  of 
an  arbitrary  point,  op2  is  the  second  Piola  Kirchhoff  stress 
tensor,  £  is  the  Green-Lagrange  strain  tensor,  and  fh  is 
the  vector  of  body  forces.  The  second  term  in  the  preceding 
equation  can  be  recognized  as  the  virtual  work  of  the  elastic 
forces,  it  can  be  rewritten  to  define  the  generalized  elastic 
forces,  that  is 


UNCLASSIFIED 


The  6th  Asian  Conference  on  Multibody  Dynamics 
August  26-30,  2012,  Shanghai,  China 


SWs=\oR:  &d\^QTs€  (18) 

V 

where  <5e  is  the  virtual  change  in  the  nodal  coordinates 
associated  with  a  particular  ANCF  finite  element  or  a  body, 
and  Qs  is  the  vector  of  the  generalized  elastic  forces.  The 
vector  of  elastic  forces  often  takes  a  fairly  complicated 
form,  especially  in  the  case  of  plasticity  formulations,  and 
is  obtained  using  numerical  integration  methods.  The 
principle  of  virtual  work  leads  to  the  following  equations  of 
motion: 

Me  +  Q-Q=0  (19) 

where  M  is  the  symmetric  mass  matrix,  and  Qe  is  the 
vector  of  body  applied  nodal  forces. 

Knowing  the  strains,  the  soil  properties,  yield  function, 
and  the  flow  rule;  the  state  of  soil  deformation  (elastic  or 
plastic)  can  be  determined.  Knowing  the  state  of 
deformation,  the  constitutive  model  appropriate  for  this 
state  can  be  used  to  determine  the  elastic  force  vector  Q  s  . 


5  COMPARISON  OF  TWO  MODELS 

Continuum-based  and  terramechanics  soil  models 
have  been  used  in  various  vehicle-terrain  interaction  studies. 
Through  these  studies  and  the  overview  presented  above, 
one  can  develop  a  comparison  of  the  two  types  of  models. 

Unlike  continuum  based  models,  terramechanics 
models  assume  a  contact  interface  appropriate  for  wheel- 
terrain  or  track-terrain  interfaces.  These  assumed  contact 
interfaces  can  take  into  consideration  contact  between  a 
flexible  wheel  and  the  terrain  through  the  use  of  concepts 
as  effective  wheel  diameter,  elliptic  soil  contact  surface,  etc. 
Continuum  models  do  not  have  this  limitation,  they  allow 
for  modeling  contact  between  arbitrary  surfaces  while 
capturing  higher  modes  of  deformation. 

As  is  evident  from  the  above  discussion,  there  exists  a 
wide  variety  of  continuum  soil  models  each  with  its  own 
merits  and  range  of  applicability.  Models  based  on  the 
concept  of  the  yield  surface  offer  a  rational  method  for 
improvement  and  development  of  new  soil  constitutive 
models.  Until  recently,  terramechanics  models'  range  of 
applicability  required  models  with  characteristic 
dimensions  similar  to  that  of  an  average  vehicle.  With 
certain  modifications  terramechanics  models  can  be  made 
relevant  to  the  simulation  of  vehicles  and  robots  of  much 
smaller  dimensions. 

While  there  are  plenty  of  differences  between 
terramechanics  and  continuum-based  soil  models,  there 


also  exist  some  similarities  between  them.  For  example, 
both  models  can  be  fitted  to  experimental  data.  Standard 
techniques  have  been  established  for  determining  the 
coefficients  that  appear  in  both  model  types.  However, 
some  continuum  models  require  extensive  testing  in  order 
to  generate  appropriate  coefficients  to  fit  experimental 
results.  Also,  both  categories  of  models  can  be  modified  to 
capture  viscous  effects. 

While  this  paper  focuses  on  discussing  the  basic 
differences  between  the  terrmechanics  and  the  continuum 
based  soil  models,  numerical  studies  that  compare  the  two 
methods  will  be  presented  by  the  authors  in  future 
investigations. 

6  SUMMARY 

In  this  paper,  several  simple  models  including 
Bekker’s  model  are  discussed.  Bekker’s  model  as  well  as 
other  parametric  and  analytical  terramechanics  models  have 
been  used  in  the  study  of  track/soil  interaction  and  can  be 
implemented  in  MBS  algorithms  using  simple  discrete 
force  elements.  More  general  continuum  plasticity  soil 
formulations  are  reviewed.  The  absolute  nodal  coordinate 
formulation  (ANCF)  which  allows  for  the  study  of 
vehicle/soil  interaction  is  also  discussed. 
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